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•  ABSTRACT 

y  A  computational! algorithm  is  developed  for  finding  the  stationary  values 
of  the  function  x' Cxj(x' Ax)^(x'Bx)^  where  A  and  B  are  positive  definite  and  C 

\  A 

is  a  symmetric  matrix.  The  square  of  the  function  under  consideration  is  the 
product  of  two  Raleigh  coefficients  x'Cx/x'Ax  and  x'Cx/x'Bx.  The  general  prob¬ 
lem  occurs  in  multivariate  analysis  in  the  computation  of  homologous  canonical 
variates  in  studying  relationships  between  two  sets  of  homologous  measurements. 
The  special  case  with  C  =  I  occurs  in  designing  control  systems  with  minimum 
norm  feedback  matrices. 
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1.  INTRODUCTION 

If  A  is  a  positive  definite  matrix  and  C  is  a  symmetric  matrix  of  order 
p,  then  it  is  well  known  that  the  stationary  values  of  the  Raleigh  coefficient 
x'Cx/x'Ax  are  the  eigen  values  of  C  with  respect  to  A  (Rao,  1973,  p.  74).  In 
particular,  if  ordered  eigen  values,  then 

X,  =  max  x'Cx  x  =  min  x'Cx  .  (1.1) 

^  X  x'Ax  *  P  X  *'Ax’ 

In  this  paper,  we  consider  the  problem  of  obtaining  the  stationary  values  of 

x'Cx  (1.2) 

(x'Ax)^(x'Bx)^ 

where  A  and  B  are  positive  definite  matrices  and  C  is  a  symmetric  matrix  of  order 
p.  The  square  of  (1.2)  is  the  product  of  the  two  Raleigh  coefficients  (x'Cx/ 
x'Ax)  and  (x'Cx/x'Bx). 

The  special  case  of  (1.2)  with  C  *  I  originally  arose  in  attempts  to  design 
control  systems  with  minimum  norm  feedback  matrices  (Kouvaritakis  and  Cameron, 
1980;  Cameron  and  Kouvaritakis,  1980)  and  also  in  the  study  of  the  stability 
of  multivariable  non-linear  feedback  systems  (Cameron,  1983).  The  general  case 
of  (1.2)  occurs  in  the  analysis  of  familial  data  when  multiple  homologous  meas¬ 
urements  are  available  on  say  father  and  son  and  the  object  is  to  determine  a 
linear  combination  of  the  measurements  which  has  the  maximum  parent-offspring 
correlation.  In  this  case,  the  dispersion  matrix  of  (Y,Z),  the  vectors  of  p 
homologous  measurements  on  father  and  son,  can  be  written  as 

A  C,  \ 

^  .  (1.3) 

C'j  B  j 

The  correlation  between  two  homologous  linear  functions  x' Y  and  x' Z  is  then 
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(1.4) 


and  the  problem  is  one  of  maximizing  or  minimizing  (1.4)  over  x<rP.  We  call 
such  optimizing  linear  functions  homologous  canonical  variates  (HCV's). 

To  obtain  the  stationary  values  of  (1.2),  we  equate  the  derivative  of  (1.2) 
with  respect  to  x  to  the  zero  vector  (Rao,  1973,  p.  72).  This  yields  the  equa¬ 
tion 

x'Cx  .  x'Cx 

7^  ^  F17  (1-5) 

which  can  be  written  in  the  equivalent  form 

XAx  +  yBx  =  2Cx 

(1.6) 

XX' Ax  =  X'  Cx  ^ 

introducing  two  additional  variables  X  and  p,  or  in  the  form 

x{A+vB)x  =  2Cx 
x'Ax  =  vx'Bx 

introducing  two  additional  variables  X  and  v. 

Since  A  and  B  are  positive  definite  matrices  there  exists  a  nonsingular 
transformation  S  such  that  A  =  SAS*  and  B  *  SS'  where  A  is  a  diagonal  matrix 
(Rao,  1973,  p.  41).  Then  writing  x  for  S' x  and  C  for  S"^C(S"^)'  ,  the  equation 
(1.7)  assumes  the  simpler  form 

x(A4vI)x  =  2Cx 
x'Ax  =  vx'x. 

If  52,...,5p  are  the  diagonal  elements  of  A  and  X2,...,Xp  are  the  components 
of  X,  then  eliminating  x  and  v  from  (1.8),  we  have  the  equations  for  X2,...,Xp 
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2x'  x  [(elCx)xj-{e^Cx)x^.  ]  =  x^x^. (s^  -6^  )x'Cx, 


(1.9) 


1  “  l»...»p» 

where  e.  is  the  elementary  vector  with  unity  as  the  i-th  component  and  zeroes 
elsewhere.  In  (1.9)  we  have  (p-1)  quartic  equations  in  (p-1)  ratios  (xg/Xj), 

. • • »(Xp_i/Xi) .  The  solution  of  these  equations  is  in  general  not  easy  except 
in  the  case  of  p  =  2  when  there  is  only  one  quartic  equation  as  observed  by 
Kouvaritakis  and  Cameron  (1980). 

In  this  paper,  we  provide  a  computational  algorithm  for  solving  the  equations 
(1.7)  in  the  general  case.  The  computer  output  gives  all  the  solutions  of  (1.7) 
and  the  corresponding  stationary  values  of  (1.2). 

2.  THE  CASE  WHERE  ALL  THE  MATRICES  ARE  DIAGONABLE 
When  all  the  matrices  A,  B  and  C  are  diagonable  by  a  common  transformation, 
the  equation  (1.8)  reduces  to 


2Fx  =  X4x  +  ux 
X'  Fx  =  XX*  Ax 


(2.1) 


where  F  is  a  diagonal  matrix  with  say  diagonal  elements.  In 

terms  of  the  components  of  x,  the  first  equation  in  (2.1)  can  be  written  as 


■jX..  -  (X6^  +  y)Xj  ,  i  -  l,...,p. 


(2.2) 


There  can  be  several  types  of  solutions  to  (2.2). 

(1)  X  =  e^.  satisfies  (2.1)  with  x  =  and  p  =  f^.  giving  the  static 

ary  value  . 

(2)  There  can  be  solutions  of  the  form  x  =  ae.  +  be..  In  such  a  case 

^  J 

2f.  =  X6^  +  p  ,  2fj  =  X6j  +  p  ,  6.  and  f.  ^  (2.3) 
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giving 


X  =  2(t,-fj)/(6^-5.),  „  =  2(f,«j-f^6,)/(S.-S,). 


(2.4) 


A  solution  of  the  form  ae.  +  be.  exists  only  if  v  =[(f  .(S  .-f  .6.)/(f  .-f.)]  £ 

I  J  I  ^  J  I  J  I 

(6^,6.).  If  this  happens,  then  x  =  ae.+be.  are  solutions  to  (2.1),  where  (a/b) 

=  [( v-'S  .)/(v-6.)]^,  yielding  the  same  stationary  value  Xv^. 

J  ^ 

(3)  There  can  be  solutions  of  the  form  ae^.  +  be^  +  but  they  lead  to 
the  same  stationary  values  as  in  (2). 

2 

An  interesting  case  is  that  of  Kantarovich  (1948),  where  6^.  =  x^. ,  f.j  =  x^. 

giving  the  stationary  values  1  corresponding  to  solutions  of  type  (1),  and  2/X^Xj/ 

(x.  +  X.)  corresponding  to  solutions  of  type  (2).  In  this  case  the  largest  value 
^  3 

is  1  and  the  smallest  is  2/nx7/(Xi  +  x„)  whero  X,  =  max{x,...,x  }  and  x  =  min 

1  p  i  p  1  p  p 

{Xj,...,Xp},  which  gives  the  celebrated  inequality  of  Kantarovich. 

Thus,  when  all  the  matrices  A,  B  and  C  are  simultaneously  diagonable,  we 

have  a  closed  form  solution  to  the  optimization  problem.  Otherwise  the  solutions 

to  (1.7)  have  to  be  obtained  through  a  suitable  algorithm  which  we  develop  in 

the  next  section. 


3.  COMPUTATIONAL  ALGORITHM  IN  THE  GENERAL  CASE 
Let  us  consider  the  basic  equation  (1.5)  in  the  form  (1.7) 


2Cx  =  x(A  +  V  B)x 


X ' Ax  =  V  X ' Bx 


(3.1) 


where  we  recall  that  A  and  B  are  positive  definite  and  C  is  a  symmetric  matrix 
all  of  order  p.  From  the  second  equation  in  (3.1),  we  find  that  v  £  [vp.vj]» 
where  and  Vp  are  the  largest  and  smallest  eigen  values  of  A  with  respect  to 
B.  For  any  given  v  £  [vpMj].  the  first  equation  in  (3.1)  provides  p  eigen  values 
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Xj{v}  ^p(v)  (3.2) 

of  2C  with  respect  to  A  +  <>6,  and  p  associated  eigen  vectors 

)»•••  »Xp (v  )  .  (3.3) 

The  pair  (v  ,x^.(v))  will  be  a  solution  of  (3.1)  if  and  only  if 

^_xi(v)A^^  (3.4) 

x!(v)Bx^.(v) 

Our  computational  algorithm  is  basically  a  search  for  v  and  a  suitable  eigen 
vector  x^. (v)  such  that  (3.4)  holds.  The  complexity  of  the  algorithm  depends 
on  the  nature  of  the  p  eigen  val ue  functions 

X..('v)»'V6  [Vp»V2]»i~l»...jP  (3.5) 

each  of  which  is  a  continuous  function  of  v  (see  Kato,  1980,  Chapter  2  for  vari¬ 
ous  results  used  in  this  section). 

If  the  rank  of  C  is  s  <  p,  then  (p-s)  functions  in  (3.5)  identically  van¬ 
ish.  All  the  solutions  of  Cx  =  0  with  X  =  0  (i.e.,  eigen  vectors  of  C  corresponding 
to  its  zero  eigen  value)  satisfy  (3.1),  and  the  stationary  value  of  (1.2)  cor¬ 
responding  to  each  such  solution  is  zero.  Then  there  is  a  fixed  number  Sj  of 
the  functions  (3.5)  such  that 

X,  (v  )  >. .  .>  X.  (v  )  >  0  (3.6) 

i  -  -  Sj 

and  a  fixed  number  S2  such  that 

X_(v)  <•••<  X_  -  .t(v  )  <  0.  (3.7) 

Let  us  start  with  X2(v).  If  Xj(v)  f  )  for  all  v  ,  then  there  is  a  unique 
continuous  eigen  vector  function  X2(v  )  associated  with  X2('v).  We  can  then  con- 
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Struct  the  continuous  function 

fl(v)  =  [x|(v)  AXj(v)/x'(v)  B  Xj(v)]  -V  (3.8) 

which  is  >  0  when  v  =  v  .  and  <  0  when  v  =  v  .so  that  there  is  at  least  one 
—  ram  —  max 

value  of  V,  say  Vj,  which  makes  (3.8)  vanish  and  provides  the  solution  [vj,Xj(vj)] 
to  (3.1).  There  may  be  more  than  one  solution  to  the  equation  fj(v)  =  0,  each 
of  which  leads  to  a  solution  of  (3.1).  Since  the  value  of  fj(v)  for  any  given 
V  is  uniquely  computable,  the  solutions  of  fj(v)  =  0  can  be  easily  found  through 
a  suitable  computer  program.  We  then  consider  )  ^nd  if  X2(v )  f  ^3^'^ ) 
any  V ,  then  the  above  procedure  can  be  implemented  leading  to  additional  solu¬ 
tions.  Now  we  go  to  additional  solutions  and  so  on.  Thus  in  the 

case  when  the  eigen  value  functions  X^.(v)  are  distinct  (no  two  meet  anywhere) 
all  the  solutions  can  be  obtained  by  considering  the  individual  ordered  eigen 
value  functions.  This  is  probably  the  case  which  often  arises  in  practice  lead¬ 
ing  to  at  least  p  solutions  of  (3.1).  Otherwise  we  proceed  as  follows. 

The  above  procedure  can  be  implemented  starting  with  Xj(v)  so  long  as  two 
successive  eigen  functions  do  not  meet.  Let  us  suppose  that  at  the  i-th  stage 
we  first  encounter  the  case 

X^(v)  =...=  )  (3.9) 

for  some  value  of  v.  Associated  with  this  repeated  root,  there  are  h  eigen  vec¬ 
tors  which  may  be  written  as  columns  of  a  matrix 

X.  =  (  .(v)  :...:  ^•+h.i(v)).  (3.10) 

Note  that  the  choice  of  the  individual  vectors  in  (3.10)  is  not  unique,  but  any 
choice  would  generate  the  same  eigen  space.  We  then  form  the  matrices 


E,  -x;  AX,,  F,  .Xj  BX, 


(3.11) 


i 
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and  find  the  largest  and  smallest  eigen  values  and  of  E.  with  respect  to 
F^.  and  the  associated  eigen  vectors  and  y^.  Then  v  will  be  a  solution  iff 
{aj-v)(a2-v)  £  0.  If  this  happens, 

V.  X  =  X.(c^yj±C2y2)  (3.12) 

2  2 

are  solutions  to  (3.1),  where  c^(a^-v)  =  C2(v-aj^),  leading  to  the  same  station¬ 
ary  value  X^.  (v).  If  (otj-v ) (a2-v )  >  0,  then  v  is  not  a  solution. 

Having  noted  the  computational  procedures  involved  in  testing  whether  a 
given  v  is  a  solution  or  not  depending  on  the  multiplicity  of  the  roots  of 

|2C-x{A+Bv)|  =0  (3.13) 

we  make  a  few  remarks  on  the  complexity  of  the  problem  one  may  run  into.  From 
the  results  in  perturbation  theory  of  symmetric  operators  it  is  known  that; 

(1)  The  number  of  distinct  roots  of  (3.12)  are  the  same  for  every  v  except 

for  a  finite  number  of  "exceptional  values"  in  where  it  can  be  less. 

(2)  The  eigen  value  functions  Xj(v)  l...iXp(v)  are  well  behaved  (holomor- 
phic)  in  the  intervals  between  the  "exceptional  points."  In  each  such  interval 
some  consecutive  eigen  functions  may  coincide,  and  the  set  of  identical  eigen 
functions  may  be  different  in  different  intervals. 

We  consider  some  examples  and  make  some  general  remarks  in  the  next  section. 


4.  ILLUSTRATIVE  EXAMPLES 


'  To  illustrate  the  computations,  we  first  consider  the  Kantarovich  problem 

where  all  the  matrices  can  be  chosen  to  be  diagonal: 


.v.\  s^V- 


V 


(4.1) 


In  this  case,  v  ^  [4,25],  where  4  is  the  smallest  and  25  is  the  largest  eigen 
value  of  A  with  respect  to  B.  The  graphs  of  X^Cv),  X2(v)  and  X2(v),  the  eigen 
value  functions  of  2C  with  respect  to  A  +  vB,  are  shown  in  Figure  1. 

[Here  Figure  1] 


We  note  that  there  are  three  possible  exceptional  points  at  which  repeated  roots 
occur  where  the  type  of  computations  described  in  (3.9)  -  (3.12)  have  to  be  done. 
Further,  between  the  exceptional  points  the  eigen  value  functions  are  distinct 
and  well  behaved. 

The  next  step  is  to  locate  the  exact  values  of  the  exceptional  points,  i.e., 
where  x^(v)  =  X2(v),  and  also  the  values  of  v  at  which  f]_(v)  vanishes  between 
the  exceptional  points.  [Note  that  fj^(v)  is  uniquely  defined  between  the  excep¬ 
tional  points  as  in  (3,8).]  This  can  be  done  by  tabulating  X2/Xj^  and  fj(v)  at 
short  intervals  of  v,  locating  the  intervals  in  which  fj^(v)  *  0  or  X2(v)/Xj^(v)  =  1, 
and  find  the  values  of  v  where  equalities  are  attained  through  a  suitable  pro¬ 
gram  for  finding  the  roots. 

If  V  is  found  such  that  f^(v)  =  0,  then  (v,  Xj(v))  is  a  solution  giving 
the  stationary  value  AT  x  (v). 

If  V  is  found  such  that  Xj(v)  =  X2{v),  then  we  have  two  eigen  vectors  say 
Xj,X2  associated  with  this  repeated  root.  We  compute  the  matrices  = 

(Xj  :  X2)'  A  (Xj  :  X2)  and  F^  =  (x^  :  X2)'  B  (x^  :  X2)  which  provide  two  eigen 
values  ^<X2  and  the  associated  eigen  vectors  ^^*^2  respect  to  Fj. 

If  (a^-v)(a2-v)  <.0,  then  v  is  a  solution  giving  the  stationary  value  AT  x^(v). 

The  vectors  at  which  this  value -is  attained  are 


V*  ’ 


1-"  IT  " 


•w’  V' 
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X  =  (x^  :  x^){c^y^±c^y^) 

2  2 

where  c^Ca^-v)  =  C2(v-a2).  If  (o2-v)(a2-v)  >  0,  then  v  is  not  a  solution. 

We  then  proceed  to  the  next  eigen  value  function  X2(v)  and  locate  the  val¬ 
ues  of  V  at  which  X2(v)/X2(v)  =  1,  and  the  non-exceptional  values  of  v  at  which 
f2(v)  =  0  and  repeat  the  above  analysis.  Finally,  we  consider  X^Cv)  and  locate 
the  non-exceptional  values  of  v  at  which  f2('^)  =  0-  The  final  tabulation  lead¬ 
ing  to  the  stationary  values  of  the  function 

(x'Cx)/(x'Ax)^(x'Bx)^  (4.2) 


is  as  follows. 


Table  1.  Stationary  values  of  (4.2) 


VJ 

x  - 

vector 

Stationary 

value 

4 

1 

0 

0 

1 

16 

0 

1 

0 

1 

25 

0 

0 

1 

1 

8* 

1//2 

1//2 

0 

.9428 

10* 

If^Z 

0 

1//2 

.9036 

20* 

0 

1//2 

1//2 

.9938 

♦Exceptional  points 


The  X  -  vectors  are  standardized  to  have  unit  nomi.  The  graphs  of  fj^(v),  f2(v) 
and  f3('')  are  shown  in  Figures  2,  3  and  4.  Note  the  discontinuities  of  each 
function  at  the  exceptional  points. 

[Here  Figures  2,  3,  and  4] 

The  next  example  is  concerned  with  the  evaluation  of  what  we  call  homolo¬ 
gous  canonical  variates  (HCV).  The  following  table  gives  the  correlation  matrix 
of  the  measurements  on  head  length  (HL),  head  width  (HW),  face  width  (FW)  and 


stature  (St)  taken  on  father  and  son.  The  problem  is  to  find  a  linear  function 
of  the  four  measurements  which  shows  the  highest  correlation  between  father  and 
son. 

Table  2.  Correlation  Matrix 


S 


.219  1.0 


0.341 

0.145 

0.243 

0.055 

0.194 

0.045 

0.066 

0.248  ' 

0.057 

-0.033 

0.111 

0.028 

0.174 

0.181 

0.187 

0.581 

<^1 

The  function  to  be  maximized  is 


p(x)  =  (x'Cx)/(x'Ax)*^(x'Bx)^S  C  =  Js(C^+C|) 


(4.3) 


In  this  case,  the  four  eigen  value  functions  are  distinct  and  corresponding  to 
each  function  there  is  only  one  root.  The  stationary  values  of  the  correlation 
function  (4.3)  and  the  standardized  vectors  at  which  they  are  attained  are  given 
in  Table  3. 
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Table  3.  Homalagou?  Canonical  Variates 


Stationary 
Value  of  p(x) 

HL 

HW 

X 

FW 

St 

.5874 

-.0076 

.1236 

-.1209 

.9826 

.3564 

.9549 

.1060 

.1679 

-.2207 

.1675 

-.1237 

-.5772 

.8001 

.1062 

-.0949 

-.3531 

.8841 

.1361 

-.2742 

5.  CONCLUDING  REMARKS 

In  practical  problems,  the  following  situations  may  arise. 

(1)  The  matrices  A,  B  and  C  are  simultaneously  diagonable  in  which  case  closed 
form  expressions  are  available  as  discussed  in  Section  2. 

(2)  The  eigen  value  functions  Xj(v ) , . . .  ,v  p(v  )  are  distinct  (no  two  have  a  com¬ 
mon  point)  in  V  £  (y  V  ]  in  which  case  the  method  described  in  the  paragraph 
containing  the  equation  (3.8)  is  applicable  leading  to  at  least  p  solutions. 

This  is  probably  the  simplest  and  the  most  frequent  case. 

(3)  The  eigen  value  functions  are  distinct  except  at  a  finite  set  of  "exceptional 
points."  In  such  a  case,  the  exceptional  points  are  dealt  with  as  in  (3.9)  - 
(3.12)  and  the  non-iHtfers^tiirg  eigen  functions  between  the  exceptional  points 

are  treated  as  in  (2)  above,  except  that  each  eigen  value  function  in  a  sub-in¬ 
terval  may  not  yield  a  root. 

(4)  A  complicated  situation  is  when  some  of  the  eigen  value  functions  coincide 
in  intervals  between  exceptional  points.  Note  that  the  number  of  distinct  eigen 
functions  in  each  interval  will  be  the  same,  although  the  eigen  functions  (v) 
that  coincide  may  be  different  in  different  intervals.  When  a  number,  say  h, 

of  eigen  value  functions  coincide  in  an  interval,  we  tabulate  aj(v)  and  otj^(v) 
defined  in  (3.11)  at  a  number  of  points  within  the  interval  and  locate  the  roots 
if  any  by  considering  the  product  (aj(v)-v)(a^(v)-v).  For  distinct  eigen  value 
functions  within  an  interval  the  procedure  indicated  in  (2)  is  followed. 


The  computational  algorithm  developed  in  this  paper  has  been  implemented 
through  FORTRAN  program  using  the  standard  routines  for  eigen  value  computations 
and  iterative  methods  for  determining  the  roots  of  an  equation  with  one  variable. 
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